The structure of coevolving infection networks 



Stefan Wieland\ Tomas Aquino^'^ and Ana Nunes^ ^'^^ 

^ Centra de Fi'sica da Materia Condensada and Departamento de Fi'sica, Faculdade de Ciencias da Universidade de Lisboa, 
P-1649-003 Lisboa, Portugal 

^ Centro de Matemdtica e Aplicagoes Fundamentais 

Faculdade de Ciencias da Universidade de Lisboa, P-1649-003 Lisboa, Portugal 



PACS 5 . 1 . Gg - Stochastic analysis methods 
PACS 8 7 . 1 . Mn - Stochastic modehng 

PACS 8 9 . 7 5 . Fb - Structures and organization in complex systems 

Abstract. - Disease awareness in infection dynamics can be modeled with adaptive contact networks whose 
rewiring rules reflect the attempt by susceptibles to avoid infectious contacts. Simulations of this type of 
models show an active phase with constant infected node density in which the interplay of disease dynamics 
and link rewiring prompts the convergence towards a well defined degree distribution, irrespective of the initial 
network topology. We develop a method to study this dynamic equilibrium and give an analytic description 
of the structure of the characteristic degree distributions and other network measures. The method applies 
to a broad class of systems and can be used to determine the steady-state topology of many other adaptive 
networks. 



Introduction. - Research on the influence of network 
topology on the global dynamics of systems composed of many 
• interacting units has been very active for the last ten years. 
\ More recently, adaptive and coevolving networks, that is, net- 
works whose topology changes along with the dynamics, have 
\ begun to be considered on two standings. The first one is to 
I explore the consequences of this interplay from the viewpoint 
■ of global dynamics f\\. The second is to understand what sort 
. of interactions may give rise to certain network architectures 
' Q. Understanding the origin of degree distributions and cor- 
: relations of real networks was first addressed by considering 
networks that evolve independently of the dynamics jS), and 
it is natural to extend this idea to networks whose evolution is 
determined by the dynamics they support. 

Infection dynamics is a natural setting to investigate the cou- 
pled evolution of the interacting elements and the network of 
interactions, because disease awareness translates into suscep- 
tibles that try to evade infection by changing their contact pat- 
terns according to the disease status of their neighbors. The 
SIS (susceptible, infected, susceptible) and the SIR (suscepti- 
ble, infected, recovered) models are the simplest models for the 
spread of epidemic infections Q. They are based on the idea 
that an infected individual infects susceptibles at a certain rate 
for a certain period, and then recovers and becomes suscepti- 
ble again (if the disease does not confer immunity), or recovers 
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and becomes immune. Stochastic versions of these models on 
networks can be mapped onto bond percolation jS)- 

Simple epidemic models like these exhibiting trivial global 
dynamics on static networks have been shown to display inter- 
esting behavior on adaptive networks, in which links between 
nodes may be created or removed according to rewiring rules 
dictated by disease awareness. In the network scheme of SIS 
dynamics, susceptible nodes {S-nodes, a fraction [S] of the to- 
tal number of nodes) are infected with rate p along links con- 
necting them to infected nodes {I-nodes, the remaining frac- 
tion [/] = 1 — [S]), whereas the latter recover with rate r. 
Gross et al. @ incorporated disease awareness into the con- 
ventional SIS model by introducing a topology-changing, yet 
link-number preserving, rewiring mechanism: S-nodes try to 
evade infection by retracting links from infected neighbors with 
rate w, and rewiring them to randomly selected S-nodes. The 
description of the model in the framework of the standard pair 
approximation yields a complex phase diagram. In addition 
to a frozen phase (where the disease dies out) and to a sta- 
tionary active phase (where [/] 7^ is constant and [AB], the 
number per node of links connecting nodes of type A and B, 
A,B<E {S,I}, is also constant), other active phases emerge 
that are absent in the SIS model without rewiring, namely an 
oscillatory phase and a bistable phase |]7]- Simulations in the 
stationary active phase show the states of the nodes and the 
links coevolving to produce and maintain a dynamic network 
topology characterized by well defined degree distributions not 
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only for the global network, but also for the subsets of S- and 
I-nodes. This remarkable fact and the structure of the charac- 
teristic degree distribution were left mostly unexplained. In the 
wake of this first paper based on SIS infection dynamics, sev- 
eral contributions have addressed extensions and modifications 
|[8l of the original model, focusing mostly on the global dynam- 
ics and the conditions for disease persistence. In a recent pa- 
per 191, steady-state degree distributions from individual -based 
simulations were shown to be very well approximated by the 
outcome of long-term numerical integration of a compartmen- 
tal model for the evolution of the node status together with the 
number of its infected and susceptible neighbors. This model 
provides an alternative to simulations to produce approximately 
the observed network topology, but it is not amenable to ana- 
lytic treatment. Overall, the description of the network topol- 
ogy underlying the steady states of these models has received 
comparatively little attention. Here we develop an analytic 
method to compute degree distributions and other measures that 
characterize this network topology. 

The node-cycle model and its solution. - We illustrate 
the method taking the original model |6| as an example of an 
adaptive network, and focusing on the stochastic process that 
each node and its links follow as the former undergoes status 
change from susceptible to infected, and gains (in the suscepti- 
ble S-stage) or loses (in the infected I-stage) links. This process 
can be treated analytically, and the stationary degree distribu- 
tions for the full node cycle (NC), from susceptible to infected 
and back again, can be obtained and compared with the degree 
distributions of the susceptible and infective subnetworks of the 
original network model. The NC approach is based on the er- 
godic properties of the network stationary state, which are a 
consequence of the random rewiring process that ensures mix- 
ing. Indeed, we observed in simulations of the steady state of 
the network model that recording ensemble statistics yields the 
same distributions as sampling a single node over a sufficiently 
long time. 

In each stage of the NC we keep track of the numbers x and 
y of susceptible and infected neighbors of the node (its joint 
degree (x,y)). Due to infection and rewiring processes, these 
numbers change at certain rates, defining a time-homogeneous 
Markov process that can be described in each stage as a finite 
random walk on a degree grid spanned by x and y (Fig. [T]i. 
The time-continuous random walks performed in each stage are 
one-step processes in both coordinates. They are coupled by 
stage transition, which takes place at constant rate r from the 
I-stage to the S-stage, and at rate py from the S-stage to the 
I-stage. The probability [x, y] = Ps{x, y, t\xo,yo) of a walker 
in the S-stage, having started at coordinates (a;o,yo), being at 
{x, y) at time t obeys the Master Equation 



I-stage 



d[x, y] 
dt 



= + r){(y + I) [x ~ l,y + I] - y[x, y]} 

" P y[x, y]+w {[x - 1, ?/] - [x, y]) 
+ Ps{{x + 1) [x + l,y- 1] - x[x,y]} 



(1) 




S-stage 



Figure 1: (Color online) Status and joint degree evolution of a node 
going through the S-stage (blue dashed line) and I-stage (red solid 
line), starting with joint degree {xo,yo) in the S-stage. Bold arrows 
show all possible transitions and their rates. 



swapping of an infected neighbor with a susceptible neighbor 
due either to rewiring or to recovery of the infected neighbor 
The second term on the right represents the transition to the I- 
stage and an overall probability mass loss as in the active phase 
all susceptible nodes eventually become infected. The stochas- 
tic dynamics of the local random variables x and y is coupled 
to the global network dynamics through the total degree gain 
rate w and through the force of infection ps that determines 
the infection rate of susceptible neighbors of the node. These 
two transitions are represented by the third and fourth terms 
on the right-hand side of Eq. ([T]i. A similar Master Equation 
holds for the I-stage, where instead of ps the rate pj accounts 
for the force of infection on susceptible neighbors of an I-node. 
We assume the correspondence parameters w, ps and pi to be 
constant and will assign values later on. 

The Master Equation ([T]i can be solved using the generat- 
ing function formalism IfTOl . Multiplying both sides of Eq.Q 
by x^7^ and summing over x and y yields for the probability 
generating function 

oo 

Fs{x^l,t\xo,yo) = ^ x'l'"Ps[x,y,t\xo,yQ) 

x,y=0 



a linear first order PDE 

dFs 



dt 



{{w + 7')x - {w + r + p)j} 



dFs 



+ Psil - x)^^ +uiix ~ 1) ■ 
OX 



(2) 



Its solution, given that Fs{Xi 7, Ojxo, yo) = X^"7^°' is 

Fsix, 1, t\xo, yo) =e-'^{^=«^+^«(*)^+^^*> 

X {ci {t)x + C2 {t)jr° 
X {C3 {t)x + Ciit)lV% 



(3) 



with the boundary conditions [— 1, y] = [x, —1] =0. The first 
term on the right hand side of Eq. ([U takes into account the 



where Cj{t),j S {1, 6} consist of linear combinations of 
exponential functions of time and cj > 0. 

The Taylor expansion of Eq. Q up to arbitrary order is easy 
to compute, yielding Ps{x,y,t\xo,yo). With the remaining 
probability mass at time t being Fs{l, 1, t\xo, yo), the normal- 
ized total occupation time at {x, y), having started at (a;o, yo), 
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reads as 



Psix,y\xo,yo) = 



J Ps{x,y,t\xo,yo)dt 


oo 

/ Fs{l,l,t\xo,yo)dt 



Given that the rate of reinfection of a S-node with y infected 
neighbours is p y, we calculate 



Ps{x,y\xo,yo) = py J Psix,y,t\xo,yo) dt 



to obtain the probability of, having started at {xo,yo), ending 
the walk in the S-stage at {x, y). In a similar way, a closed- 
form expression for the I-stage can be obtained, where Pi = Pi 
because of recovery being coordinate-independent. 

The elements of Pa, A e {S, /}, are the transition probabil- 
ities that map a random walker's initial coordinate distribution 
(ICD) in stage A, $^(a;, y), onto the ICD of the subsequent 
stage according to 

CJO 

'^i,six,y)= ^ Ps,iix,y\xo,yo)'^s,i{xo,yo)- (4) 

xo,yo=o 

The ICD of the S-stage after a full node cycle changes accord- 
ing to 

oo 

*!s(2;,y)= X! Pi^>y\^o,yo)'^si.xo,yo), 

xo,yo=o 

where the P{x,y\xo,yo) are given by the Chapman- 
Kolmogorov identity 



P{x,y\xo,yo) 



E 

x' ,y'—0 



Pi{x, y\x', y')Ps{x', y'\xo, yo) 



Setting a cutoff for the maximum degree, the elements of P en- 
code a finite ergodic Markov chain, and therefore an arbitrary 
ICD in the S-stage will converge upon iteration of the node cy- 
cle to a unique stationary state $g (x, y) which, according to the 
Perron-Frobenius theorem, can be found as the right eigenvec- 
tor of P associated with the eigenvalue 1. Once $J j{x, y) are 
known, the normalized total occupation times of coordinates in 
the stationary state can be computed in either stage as 

oo 

Ps,ii^^y)= E Psjix>y\^o,yo)'^*sjixo,yo), (5) 



where in particular 



P;(x,2/) = $*(a:,y) 



because of Eq. (HJi and Pi ~ Pi. In each stage A, A ^ {S, I}, 
the average number of infected and susceptible neighbors, (xa) 
and (yA), and the average degree, (fc/i) = (xa) + {yA), can be 
computed from the distributions Pg j{x,y). 



The lifetime distributions in each stage Tsj{t, x, y) of ran- 
dom walkers with initial coordinates (a;, y) and the generating 
functions Fsj are related through 



FsA^AAx,y)^l~ / Tsj{t',x,y)At' , 



with overall lifetime distributions Ts.i{t) 

oo 

TsAt)= E nA^.y)Ts,i{i^^^y) 

x,y—0 

and survival functions Lsj{t) in the steady state 

t 



oo 

= E '^*sAx,y)FsA^^Ax.y)- (6) 

x,y—0 

The average duration T5 / of each stage in the stationary state 
is thus given by 

00 °? 

Ts,i = E ^*sAx,y) j Fsj{i,iAx,y)'^t, 

x,y=0 Q 

It then follows from the solution of the analogue of Eq. for 
the I-stage that t/ — l/r as expected, because there random 
walkers switch to the S-stage at a coordinate-independent con- 
stant rate r. 

Correspondence with the network model and compari- 
son with simulations. - The stationary state of the NC pro- 
cess was solved analytically for any choice of the parameters 
w, r, p and w, ps, pi- These solutions encompass a much 
larger set of systems than the original network model of |l6l, 
because the correspondence parameters, tying network to NC 
dynamics, have yet to be determined. A first estimate for these 
parameters may be obtained from the equilibrium link num- 
ber per node predicted by the pair approximation model of ||6l . 
Combining this approximation with the NC description yields 
w = w — [/]), and similar expressions hold for the other 

correspondence parameters. The results for the subensemble 
degree distributions obtained through this approach reproduce 
the qualitative behaviour of the network simulations but the 
quantitative agreement is poor (results not shown). An alterna- 
tive way of constraining the correspondence parameters comes 
out of embedding the NC description in a given network of 
fixed degree. On one hand, the correspondence parameters 
must be such that link depletion in the I-stage and link addi- 
tion in the S-stage balance out for that degree. On the other 
hand, they should ensure that computing the numbers of links 
between nodes of type S and / does not depend on whether one 
sums over the infected neighbors of all the S-nodes or over the 
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susceptible neighbors of all I-nodes. These consistency condi- 
tions will provide a self-contained determination of the corre- 
spondence parameters within the NC framework which is inde- 
pendent of the pair approximation and does not involve network 
simulations for parameter fitting. 

For the NC process to represent the stationary dynamics of 
the network, so that the Pg j{x, y) correspond to the network 
subensemble degree distributions and ts and r/ are propor- 
tional to the fractions [/] and [S] of infected and susceptible 
nodes in the network, we shall therefore impose the following 
two conditions: 

(i) The average node degree must be equal to the average degree 
ik) that is chosen a priori in the network model, 

Ts{ks) + Ti{kj) ^ {k){Ts +Ti), 

where (ks.i) are the average degrees computed from Eq. (|5]l. 

(ii) The S- and I-stages represent nodes that are each other's 
neighbors, and therefore ts {ys) = tj (xj) for the mean number 
of susceptible and infected neighbors in the two stages, with the 
averages again computed from Eq. Q. 

Computing the NC stationary state and choosing values of 
the correspondence parameters such that the consistency con- 
ditions (i) and (ii) hold, the NC description can be mapped onto 
the stationary state of a particular network model. However, the 
constraints set by (i) and (ii) may be impossible to satisfy ex- 
actly. Because of the mean field approximation involved in the 
assumption of constant ps and pj, the stochastic process we 
have dealt with does not capture the weak correlations present 
in the network. This translates into (i) and (ii) not being exactly 
fulfilled, and so the description of the stationary state given by 
the NC is only approximate. In order to compare the analytic 
results with simulations on networks, we have used the free pa- 
rameters w, Ps and pj to minimize the squared distance of the 
NC output to the target values set by conditions (i) and (ii), 
bringing the stochastic model to represent as accurately as pos- 
sible the state and degree evolution of a network node. We shall 
refer to the analytic stationary degree distributions obtained 
through the method outlined above for the optimal choice of w, 
PS and pj as the NC degree distributions. Results for a high and 
a low rewiring regime in a network of average degree (k) = 7 
are presented in Fig. |2l For different choices of w, r and p 
in the stationary active phase, the NC degree distributions as 
well as ^j{x,y) show very good quantitative agreement with 
those obtained from long Monte Carlo (MC) simulations of the 
network model of ||6l and ||9l, as do cumulative node lifetime 
distributions as given by Eq. ^ (Fig. O. Our method predicts, 
through combinatorial handling of Pg j {x, y), steady-state den- 
sities of any star motif. With steady-state node densities being 
given through 

ts/ti = (1 - [I]nc)I[I]nc. 

the density of triplets with a central I-node connected to a S- 
and another I-node, normalized by the system size, is illustra- 
tively computed as 

oo 




Figure 2: (Color online) Subensemble steady-state degree distribu- 
tions for susceptibles (circles), infected (diamonds) as well as the 
steady-state ICD for the I-stage (squares) for two different rewiring 
regimes. Insets: Linear plots of the distributions in their dominat- 
ing degree range; comparison of steady-state prevalence [/]mc taken 
from MC simulations and [I]nc computed by the NC framework. 
Solid lines are predictions by the node cycle. Parameters p — 0.008, 
r = 0.005. (a): w = 0.025, w = 0.12, ps = 0.044, pi = 0.049. (b): 
w = 0.050, w = 0.22, ps = 0.042, pi = 0.045. Cutoff for overall 
degree in NC matrices is fcmax = 80. MC simulations according to 
II II with N — 5 ■ 10* nodes, (k) = 7, initial Erdos-Renyi graph and 
initial prevalence [/]o — 0.6. Statistics were recorded att — 3- 
for 10'^ network realizations. Error bars are smaller than markers. 

For the parameters of Fig.|2a), this yields [J/S'Jnc = 3.844, 
with [IIS]mc = 3.824 ± 0.004 from MC simulations. The 
small discrepancies (see Fig. [2]i between steady-state moment 
densities observed in MC simulations and those predicted by 
the NC method are another consequence of the constant ps and 
constant pj approximation. An improvement on this approxi- 
mation would be to consider ps and pj to be linear functions 
of time chosen to represent status and degree correlations of 
the network model. For instance, pi should increase in time 
along the I-stage, because older (lower degree) I-nodes have 
susceptible neighbors that are younger than average S-nodes, 
and younger S-nodes in turn have more than the average num- 
ber of infected neighbors. Although it would no longer be pos- 
sible in this case to find a simple solution like Eq. (O for the 
generating functions Fsj, the whole NC procedure could in 
principle be carried out numerically for any set of positive rate 
functions to yield an irreducible transition matrix P, a unique 
stationary degree distribution Pg j{x, y) for each stage, and im- 
proved approximations for all the moment densities that char- 
acterize the network's steady state. 
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Figure 3: Survival function of S-nodes in steady state. Inset: Linear 

plot of the distribution for most prevalent lifetimes 

; comparison of steady-state mean lifetimes ts taken from MC 

simulations and computed by the NC framework. The solid 
line is the prediction by the node cycle. Parameters and initial 
conditions as in Fig. 2(b). Error bars are smaller than markers. 



Discussion and conclusions. - Apart from giving a com- 
plete description of the network's steady state that is in excel- 
lent quantitative agreement with the results of MC simulations, 
the NC method also shows that convergence of an arbitrary ICD 
to a unique steady-state topology determined by the parameters 
It;, r, p and by the network's fixed degree will take place once 
the rate functions associated with the correspondence param- 
eters w, ps and pj are defined. These may be simply cho- 
sen as the constants that optimize the consistency conditions 
(i) and (ii), or, as mentioned above, may be taken as functions 
of time with more free parameters, allowing a closer match of 
the NC model with the network-based system. In either case, 
P encodes a time-homogeneous Markov chain, meaning that 
the method applies only to the network model's stationary state 
and not to its transients. However, the correspondence parame- 
ters embody the mean-field approximations involved in the NC 
construction, and in the simplest case they are related only with 
two steady-state moment densities, [SI] and [/] . If, as predicted 
in 161 and confirmed by simulations, in the stationary active 
phase these densities reach equilibrium independently of net- 
work structure, then the NC approach shows that this implies 
that the network topology, as described by various probability 
distributions, wiU reach equilibrium as well. 

Hence the pairwise model of ||6| and our NC approach should 
be considered complementary analytic frameworks to describe 
approximately adaptive networks. The former allows for a 
complete global treatment of system dynamics, at the cost of 
confining the level of description to low-order moment densi- 
ties. The latter is a local model, in that it identifies active steady 
states and gives a comprehensive description of their equilib- 
rium dynamics through various characteristic distributions. 



In conclusion, taking the SIS network model with rewiring 
of Gross et al. as an example, we developed a self-sufficient 
method to derive analytically the steady-state subensemble de- 
gree distributions, and other measures that characterize the 
unique steady state in the simple endemic phase of the model. 
The only requirements of the construction are that we have a 
cyclic two-state system, that the state changes are due to node 
or contact processes with constant rates, and that the network 
dynamics is due to link processes also with constant rates. The 
requirement of constant rates may be relaxed at the expense 
of an analytic solution for the probability generating functions 
not being available in this case. As for the improved approxi- 
mations briefly discussed in the previous section, dealing with 
non-linear transition rates would require the whole NC proce- 
dure to be carried out numerically to yield the stationary de- 
gree distribution for each stage. The requirement of a two-state 
system may also be relaxed to include cyclic multi-state sys- 
tems such as SIRS or rock-paper-scissors dynamics ITSl . The 
method developed here can therefore serve as a tool to explore 
the relation between the node's status change and rewiring rules 
and the resulting network structure in different contexts, rang- 
ing from infection to opinion dynamics. 
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